## Replication Files for
##
## From Delegation to Defiance: Voter Support for Party Differentiation in Coalition Governments.
## Forthcoming in Party Politics
##
## Pit Rieger - Oct 10 2025
library(here)
library(haven)
library(tidyverse)
library(modelsummary)
library(tinytable)
library(stargazer)
library(marginaleffects)
library(fixest)
library(rdrobust)
library(egg)
library(cowplot)
library(gridExtra)
library(WeightIt)
library(cobalt)
library(dotwhisker)
library(ggpubr)
library(grid)

# utility
source(here("03_util.R"))
treatdate = as.Date("2024-06-17")

# Read data ----
  # Main
  d_combined = readRDS(here("data_combined.RDS"))
  d_combined$running_visual = d_combined$running - 0.5
  d_combined = d_combined %>%
    mutate(lr_tri_label = case_when(lr_tri == "left" ~ "Left-leaning voters",
                                    lr_tri == "center" ~ "Centrist voters",
                                    lr_tri == "right" ~ "Right-leaning voters",
                                    T ~ NA),
           lr_tri_label = factor(lr_tri_label, levels = c("Left-leaning voters", "Centrist voters", "Right-leaning voters")))

  # AUTNES all waves
  d_autnes_all = readRDS(here("data_autnes_all.RDS"))

  # Google trends
  google = read_csv(here("data_googletrend.csv"), skip = 2)
  google = google %>%
    rename(date = Tag,
           Gewessler = `Gewessler: (Österreich)`,
           Renaturierungsgesetz = `renaturierungsgesetz: (Österreich)`,
           Nehammer = `Nehammer: (Österreich)`) %>%
    pivot_longer(2:4, names_to = "term",
                 values_to = "index")

  google_long = read_csv(here("data_googletrend_long.csv"), skip = 2)
  google_long = google_long %>%
    rename(date = Woche,
           Gewessler = `Gewessler: (Österreich)`,
           Renaturierungsgesetz = `Renaturierungsgesetz: (Österreich)`,
           Nehammer = `Nehammer: (Österreich)`) %>%
    pivot_longer(2:4, names_to = "term",
                 values_to = "index")

  google_hourly = read_csv(here("data_googletrend_hourly2.csv"), skip = 2)
  google_hourly = google_hourly %>%
    rename(datetime = Zeit,
           Gewessler = `Gewessler: (Österreich)`,
           Renaturierungsgesetz = `Renaturierungsgesetz: (Österreich)`,
           Nehammer = `Nehammer: (Österreich)`) %>%
    pivot_longer(2:4, names_to = "term",
                 values_to = "index")

# Main results ----
  ## RDD plot (Fig 1) ----
  d_combined_plot = d_combined %>%
    group_by(lr_tri_label) %>%
    mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T),
           ptv_OVP = ptv_OVP - mean(ptv_OVP, na.rm = T))

  grid.arrange(RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "lr_tri_label") +
                 scale_y_continuous(limits = c(-2.4, 1.5)) +
                 labs(y = "Propensity to vote for Greens (centered)",
                      x = " "),
               RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "lr_tri_label") +
                 scale_y_continuous(limits = c(-2.4, 1.5))+
                 labs(y = "Propensity to vote for ÖVP (centered)"))

    ## RDD table (Tab 1) ----
    summary(fit1.main <- with(d_combined %>% filter(lr_tri == "left"),
                              rdrobust(ptv_GRUN, running, all = T)))
    summary(fit2.main <- with(d_combined %>% filter(lr_tri == "center"),
                              rdrobust(ptv_GRUN, running, all = T)))
    summary(fit3.main <- with(d_combined %>% filter(lr_tri == "right"),
                              rdrobust(ptv_GRUN, running, all = T)))

    summary(fit1.main_OVP <- with(d_combined %>% filter(lr_tri == "left"),
                                  rdrobust(ptv_OVP, running, all = T)))
    summary(fit2.main_OVP <- with(d_combined %>% filter(lr_tri == "center"),
                                  rdrobust(ptv_OVP, running, all = T)))
    summary(fit3.main_OVP <- with(d_combined %>% filter(lr_tri == "right"),
                                  rdrobust(ptv_OVP, running, all = T)))
    models_main = list(fit1.main, fit2.main, fit3.main,
                       fit1.main_OVP, fit2.main_OVP, fit3.main_OVP)

    FE_add = data.frame(matrix(c("Sample", "Left", "Center", "Right", "Left", "Center", "Right"),
                               ncol = length(models_main) + 1, byrow = T))
    attr(FE_add, "position") = c(7, 8)

    regtab = modelsummary(models_main,
                          output = "latex",
                          add_rows = FE_add,
                          escape=FALSE,
                          stars = T,
                          title = "RD results. Standard errors in parentheses. Dependent variables: Propensity to vote for Greens (models 1-3) and ÖVP (models 4-6). \\label{tab:rdtab_ptvGRUNOVP}")# c('*' = .05, '**' = .01, "***" = .001))
    regtab

# Appendix ----
  ## Google Trends (Fig A1) ----
    vote_week = as.Date("2024-06-16") # SIC! because the long google trends is at weekly level starting on Sunday.
    p0 = ggplot(google_long,
                aes(x = date, y = index, linetype = term)) +
      annotate(geom = "segment",
               x = vote_week,
               xend = vote_week,
               y = 0,
               yend = 110,
               color = "grey60",
               linewidth = 0.4) +
      annotate(geom = "label",
               x = vote_week, y = 110,
               label = "Council vote") +
      geom_line() +
      labs(y = "Google search frequency") +
      scale_y_continuous(breaks = c(0, 100), labels = c("Min", "Max")) +
      scale_x_date(expand = expansion(add = 0),
                   #breaks = a,
                   date_minor_breaks = "months",
                   date_labels = "%d %b") +
      guides(linetype = guide_legend(title = "Search term")) +
      theme_bw() +
      theme(axis.title.x = element_blank(),
            #        legend.key.width = unit(1/2, "cm"),
            legend.title = element_blank(),
            legend.position = "bottom",
            plot.margin = margin(6,18,6,6))
    p0

    vote = as.Date("2024-06-17")
    p1 = ggplot(google %>% filter(date >= min(d_combined$date) - 1 &
                                    date <= max(d_combined$date) + 1) ,
                aes(x = date, y = index, linetype = term)) +
      annotate(geom = "segment",
               x = vote,
               xend = vote,
               y = 0,
               yend = 110,
               color = "grey60",
               linewidth = 0.4) +
      annotate(geom = "label",
               x = vote, y = 110,
               label = "Council vote") +
      geom_line() +
      labs(y = "Google search frequency") +
      scale_y_continuous(breaks = c(0, 100), labels = c("Min", "Max")) +
      scale_x_date(expand = expansion(add = 0),
                   breaks = treatdate + seq(-6, 6, 3),
                   date_minor_breaks = "day",
                   date_labels = "%d %B") +
      guides(linetype = guide_legend(title = "Search term")) +
      theme_bw() +
      theme(axis.title.x = element_blank(),
            #        legend.key.width = unit(1/2, "cm"),
            legend.title = element_blank(),
            legend.position = "bottom",
            plot.margin = margin(6,18,6,6))
    p1
    pc = as.POSIXct("2024-06-16 14:00", tz = "GMT")
    vote = as.POSIXct("2024-06-17 10:00", tz = "GMT")

    p2 = ggplot(google_hourly %>% filter(datetime < as.POSIXct("2024-06-18", tz = "GMT")&
                                           datetime > as.POSIXct("2024-06-16 00:00", tz = "GMT")),
                aes(x = datetime, y = index, linetype = term)) +
      geom_segment(mapping = aes(x = pc,
                                 xend = pc,
                                 y = 0, yend = 110),
                   color = "grey60",
                   linewidth = 0.4) +
      geom_label(mapping = aes(x = pc, y = 110,
                               label = "Press conference")) +
      geom_segment(mapping = aes(x = vote,
                                 xend = vote,
                                 y = 0, yend = 110),
                   color = "grey60",
                   linewidth = 0.4) +
      geom_line() +
      geom_label(mapping = aes(x = vote, y = 110,
                               label = "Council vote")) +
      labs(y = "Google search frequency") +
      scale_y_continuous(breaks = c(0, 100), labels = c("Min", "Max")) +
      guides(linetype = guide_legend(title = "Search term")) +
      theme_bw() +
      theme(axis.title.x = element_blank(),
            legend.title = element_blank(),
            legend.position = "bottom",
            legend.key.width = unit(1, "cm"),
            legend.title.position = "top",
            plot.margin = margin(6,18,6,6))

    # combine plots
    mylegend <- g_legend(p0)
    p = grid.arrange(ggarrange(p0 + rremove("legend") + rremove("ylab"),
                               p1 + rremove("legend") + rremove("ylab"),
                               p2 + rremove("legend") + rremove("ylab"),
                               ncol = 1),
                     arrangeGrob(mylegend),
                     ncol = 1,
                     heights = c(8, 0.5))
    annotate_figure(grid.arrange(p),
                    left = textGrob("Google search frequency", rot = 90, vjust = 1, gp = gpar(cex = 1)))

  ## Survey information (Fig A3) ----
    d_combined_sum = d_combined %>% group_by(date, survey) %>%
      summarize(n = n())
    d_combined_sum$label = NA
    d_combined_sum$label[d_combined_sum$date == treatdate] = "Gewessler votes for\n nature restoration law"
    d_combined_sum$Survey = toupper(d_combined_sum$survey)

   ggplot(d_combined_sum, aes(x = date, y = n, fill = Survey, label = label)) +
      geom_bar(stat = "identity") +
      geom_text(inherit.aes = FALSE,
                data = d_combined_sum %>% filter(!is.na(label)) %>%
                  group_by(date) %>%
                  summarize(n = sum(n), label = first(label)),
                mapping = aes(x = date + 2, y = n + 600, label = label),  # Adjust vertical position
                vjust = 0.4,
                hjust = 0) +  # Align text above
      geom_curve(inherit.aes = FALSE,
                 data = d_combined_sum %>% filter(!is.na(label)) %>%
                   group_by(date) %>%
                   summarize(n = sum(n)),
                 aes(x = date + 2, xend = date, y = n + 600, yend = n),
                 arrow = arrow(length = unit(0.1, "inches"), type = "closed"),
                 curvature = .3,  # Adjust curvature (-ve for downward arch)
                 linewidth = 0.6,
                 color = "black") +
      scale_fill_manual(values = c("#2E86C1", "#A93226")) +
      scale_x_date(expand = expansion(add = 0.55),
                   breaks = treatdate + seq(-6, 6, 3),
                   date_minor_breaks = "day",
                   date_labels = "%d %B") +
      theme_bw() +
      theme(legend.title = element_blank(),
            axis.title.x = element_blank(),
            legend.position = "bottom") +
      labs(y = "Number of respondents")

    # field dates
    d_combined %>% group_by(survey) %>%
      reframe(range(date),
              c("min", "max"))

  ## PTV other parties (Fig A4) ----
    ### SPO
    d_combined_plot = d_combined %>%
      group_by(lr_tri_label) %>%
      mutate(ptv_SPO = ptv_SPO - mean(ptv_SPO,na.rm = T))
    RDDplot(d_combined_plot, "ptv_SPO", "running", "treat", "lr_tri_label") +
      labs(y = "Propensity to vote for SPÖ (centered)",
           x = "Days since Council vote")

    ### FPO
    d_combined_plot = d_combined %>%
      group_by(lr_tri_label) %>%
      mutate(ptv_FPO = ptv_FPO - mean(ptv_FPO,na.rm = T))
    RDDplot(d_combined_plot, "ptv_FPO", "running", "treat", "lr_tri_label") +
      labs(y = "Propensity to vote for FPÖ (centered)",
           x = "Days since Council vote")

    ### NEOS
    d_combined_plot = d_combined %>%
      group_by(lr_tri_label) %>%
      mutate(ptv_NEOS = ptv_NEOS - mean(ptv_NEOS,na.rm = T))
    RDDplot(d_combined_plot, "ptv_NEOS", "running", "treat", "lr_tri_label") +
      labs(y = "Propensity to vote for NEOS (centered)",
           x = "Days since Council vote")

  ## Government satisfaction ----
    ### RDD plot (Fig A5) ----
    d_combined_plot = d_combined %>%
      group_by(lr_tri_label) %>%
      mutate(gov_satis = gov_satis - mean(gov_satis,na.rm = T))
    RDDplot(d_combined, "gov_satis", "running", "treat", "lr_tri_label") +
      labs(y = "Government satisfaction",
           x = "Days since Council vote")

    ### RDD table (Tab A2) ----
    summary(fit1 <- with(d_combined %>% filter(lr_tri == "left"),
                         rdrobust(gov_satis, running, all = T)))
    summary(fit2 <- with(d_combined %>% filter(lr_tri == "center"),
                         rdrobust(gov_satis, running, all = T)))
    summary(fit3 <- with(d_combined %>% filter(lr_tri == "right"),
                         rdrobust(gov_satis, running, all = T)))
    models_main = list(fit1, fit2, fit3)

    FE_add = data.frame(matrix(c("Sample", "Left", "Center", "Right"),
                               ncol = length(models_main) + 1, byrow = T))
    attr(FE_add, "position") = c(7, 8)

    # modelsummary
    regtab = modelsummary(models_main,
                          output = "latex",
                          add_rows = FE_add,
                          escape=FALSE,
                          stars = T,
                          title = "RD results. Standard errors in parentheses. Dependent variable: Government satisfaction.\\label{tab:rdtab_govsatis}")# c('*' = .05, '**' = .01, "***" = .001))
    regtab

  ## Political Interest ----
    ### Greens (Fig A6) ----
    d_combined_plot = d_combined %>%
      filter(polinterest == "high") %>%
      group_by(lr_tri) %>%
      mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T))
    p1 = RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "lr_tri_label") +
      theme(plot.title = element_text(face = "bold")) +
      labs(y = "Propensity to vote for Greens (centered)",
           x = "Days since Council vote",
           title = "High political interest voters")

    d_combined_plot = d_combined %>%
      filter(polinterest == "low") %>%
      group_by(lr_tri) %>%
      mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T))
    p2 = RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "lr_tri_label") +
      scale_y_continuous(breaks = c(-2.0, -1.0, 0, 1.0, 2.0)) +
      theme(plot.title = element_text(face = "bold")) +
      labs(y = "Propensity to vote for Greens (centered)",
           x = "Days since Council vote",
           title = "Low political interest voters")
    p3 = ggarrange(p1 + rremove("ylab"),
                   p2 + rremove("ylab"),
                   ncol = 1)
    annotate_figure(p3,
                    left = textGrob("Propensity to vote for Greens (centered)", rot = 90, vjust = 1, gp = gpar(cex = 1)))


    ### OVP (Fig A7)----
    d_combined_plot = d_combined %>%
      filter(polinterest == "high") %>%
      group_by(lr_tri) %>%
      mutate(ptv_OVP = ptv_OVP - mean(ptv_OVP,na.rm = T))
    p1 = RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "lr_tri_label") +
      theme(plot.title = element_text(face = "bold")) +
      labs(y = "Propensity to vote for ÖVP (centered)",
           x = "Days since Council vote",
           title = "High political interest voters")

    d_combined_plot = d_combined %>%
      filter(polinterest == "low") %>%
      group_by(lr_tri) %>%
      mutate(ptv_OVP = ptv_OVP - mean(ptv_OVP,na.rm = T))
    p2 = RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "lr_tri_label") +
      scale_y_continuous(breaks = c(-2.0, -1.0, 0, 1.0, 2.0)) +
      theme(plot.title = element_text(face = "bold")) +
      labs(y = "Propensity to vote for ÖVP (centered)",
           x = "Days since Council vote",
           title = "Low political interest voters")
    p3 = ggarrange(p1 + rremove("ylab"),
                   p2 + rremove("ylab"),
                   ncol = 1)

    annotate_figure(p3,
                    left = textGrob("Propensity to vote for ÖVP (centered)", rot = 90, vjust = 1, gp = gpar(cex = 1)))

  ## Covariate Tests (Fig A8) ----
    covariate_ests = data.frame()

    get_estimates = function(data, y, running = "running"){
      est = rdrobust(data[,y], data[,running])
      est90 = rdrobust(data[,y], data[,running], level = 90)
      data.frame(est = est$Estimate[1],
                 lower = est$ci[1,1],
                 upper = est$ci[1,2],
                 lower90 = est90$ci[1,1],
                 upper90 = est90$ci[1,2])
    }
    get_estimates_tri = function(data, y, running = "running"){
      bind_rows(get_estimates(data %>% filter(lr_tri == "left" & !is.na(sym(y))),
                              y = y, running = running) %>%
                  mutate(lr_tri = "left"),
                get_estimates(data %>% filter(lr_tri == "center" & !is.na(sym(y))),
                              y = y, running = running) %>%
                  mutate(lr_tri = "center"),
                get_estimates(data %>% filter(lr_tri == "right" & !is.na(sym(y))),
                              y = y, running = running) %>%
                  mutate(lr_tri = "right"))
    }

    # estimate age
    covariate_ests = bind_rows(covariate_ests,
                               get_estimates_tri(d_combined, y = "age_scaled") %>%
                                 mutate(variable = "Age"))

    # estimate gender
    covariate_ests = bind_rows(covariate_ests,
                               get_estimates_tri(d_combined, y = "gender_num_scaled") %>%
                                 mutate(variable = "Gender"))

    # estimate edu
    covariate_ests = bind_rows(covariate_ests,
                               get_estimates_tri(d_combined, y = "edu2_scaled") %>%
                                 mutate(variable = "Education"))

    # estimate rural
    covariate_ests = bind_rows(covariate_ests,
                               get_estimates_tri(d_combined, y = "rural_scaled") %>%
                                 mutate(variable = "Rural"))


    covariate_ests$lr_tri_label = case_when(covariate_ests$lr_tri == "left" ~ "Left-leaning voters",
                                            covariate_ests$lr_tri == "center" ~ "Centrist voters",
                                            covariate_ests$lr_tri == "right" ~ "Right-leaning voters")
    covariate_ests$lr_tri_label = factor(covariate_ests$lr_tri_label, levels = c("Left-leaning voters", "Centrist voters", "Right-leaning voters"))
    covariate_ests$variable = factor(covariate_ests$variable,
                                     levels = rev(c("Age", "Gender", "Education",
                                                    "Rural")))

    ggplot(covariate_ests, aes(x = est, xmin = lower, xmax = upper, y = variable,
                               color = lr_tri_label, shape = lr_tri_label)) +
      geom_vline(xintercept = 0, linetype = 2) +
      geom_pointrange(position = position_dodge(width = 0.6),
                      size = 1,
                      linewidth = 1) +
      geom_linerange(inherit.aes = F,
                     mapping = aes(x = est, xmin = lower90, xmax = upper90,
                                   y = variable, color = lr_tri_label),
                     linewidth = 2,

                     position = position_dodge(width = 0.6)) +

      scale_color_manual(values = c("#2E86C1", "#A93226", "orange")) +
      theme_bw() +
      theme(legend.title = element_blank(),
            legend.position = "bottom",
            legend.key.width = unit(2, "cm"),
            axis.title.y = element_blank()) +
      labs(x = "RD estimate")

  ## Covariate adjustment ----
    ### Greens (Tab A3) ----
      # matching
      d_combined$temp_id = 1:nrow(d_combined)
      d_entropy = d_combined %>%
        select(temp_id, treat, running, gender, age_scaled, edu2, rural) %>%
        na.omit()
      fit.match = weightit(treat ~
                             gender +
                             age_scaled +
                             as.factor(edu2) +
                             as.factor(rural),
                           d_entropy,
                           method = "ebal")
      love.plot(fit.match, stars = "std")
      d_entropy$w = fit.match$weights
      hist(d_entropy$w)
      d_combined$w = NULL
      d_combined = left_join(d_combined, d_entropy %>% select(temp_id, w))

      # matching models
      fit.RDmatch1 = with(d_combined %>% filter(!is.na(w) & lr_tri == "left"),
                          rdrobust(ptv_GRUN, running, weights = w,
                                   all = T))
      summary(fit.RDmatch1)
      fit.RDmatch2 = with(d_combined %>% filter(!is.na(w) & lr_tri == "center"),
                          rdrobust(ptv_GRUN, running, weights = w,
                                   all = T))
      summary(fit.RDmatch2)
      fit.RDmatch3 = with(d_combined %>% filter(!is.na(w) & lr_tri == "right"),
                          rdrobust(ptv_GRUN, running, weights = w,
                                   all = T))
      summary(fit.RDmatch3)

      # control variables
      fit.covariate1 = with(d_combined %>% filter(!is.na(w) & lr_tri == "left"),
                            rdrobust(ptv_GRUN, running,
                                     covs = model.matrix(~ -1 +
                                                           gender +
                                                           age_scaled +
                                                           as.factor(edu2) +
                                                           as.factor(rural)),
                                     all = T))
      summary(fit.covariate1)
      fit.covariate2 = with(d_combined %>% filter(!is.na(w) & lr_tri == "center"),
                            rdrobust(ptv_GRUN, running,
                                     covs = model.matrix(~ -1 +
                                                           gender +
                                                           age_scaled +
                                                           as.factor(edu2) +
                                                           as.factor(rural)),
                                     all = T))
      summary(fit.covariate2)
      fit.covariate3 = with(d_combined %>% filter(!is.na(w) & lr_tri == "right"),
                            rdrobust(ptv_GRUN, running,
                                     covs = model.matrix(~ -1 +
                                                           gender +
                                                           age_scaled +
                                                           as.factor(edu2) +
                                                           as.factor(rural)),
                                     all = T))
      summary(fit.covariate3)

      # RDD table
      models_main = list(fit.RDmatch1, fit.covariate1,
                         fit.RDmatch2, fit.covariate2,
                         fit.RDmatch3, fit.covariate3)

      FE_add = data.frame(matrix(c("Adjustment",  "Matching", "Controls", "Matching", "Controls", "Matching", "Controls",
                                   "Sample", "Left", "Left", "Center", "Center", "Right", "Right"),
                                 ncol = length(models_main) + 1, byrow = T))
      attr(FE_add, "position") = c(7, 8)

      regtab = modelsummary(models_main,
                            output = "latex",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "RD results with covariate adjustment. Standard errors in parentheses. Dependent variable: Propensity to vote for Greens \\label{tab:rdtab_ptvGRUN_covariates.tex}")
      regtab

    ### OVP (Tab A4) ----
      # matching models
      fit.RDmatch1_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "left"),
                              rdrobust(ptv_OVP, running, weights = w,
                                       all = T))
      summary(fit.RDmatch1_OVP)
      fit.RDmatch2_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "center"),
                              rdrobust(ptv_OVP, running, weights = w,
                                       all = T))
      summary(fit.RDmatch2_OVP)
      fit.RDmatch3_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "right"),
                              rdrobust(ptv_OVP, running, weights = w,
                                       all = T))
      summary(fit.RDmatch3_OVP)

      # control variables
      fit.covariate1_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "left"),
                                rdrobust(ptv_OVP, running,
                                         covs = model.matrix(~ -1 +
                                                               gender +
                                                               age_scaled +
                                                               as.factor(edu2) +
                                                               as.factor(rural)),
                                         all = T))
      summary(fit.covariate1_OVP)
      fit.covariate2_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "center"),
                                rdrobust(ptv_OVP, running,
                                         covs = model.matrix(~ -1 +
                                                               gender +
                                                               age_scaled +
                                                               as.factor(edu2) +
                                                               as.factor(rural)),
                                         all = T))
      summary(fit.covariate2_OVP)
      fit.covariate3_OVP = with(d_combined %>% filter(!is.na(w) & lr_tri == "right"),
                                rdrobust(ptv_OVP, running,
                                         covs = model.matrix(~ -1 +
                                                               gender +
                                                               age_scaled +
                                                               as.factor(edu2) +
                                                               as.factor(rural)),
                                         all = T))
      summary(fit.covariate3_OVP)

      # RDD table
      models_main = list(fit.RDmatch1_OVP, fit.covariate1_OVP,
                         fit.RDmatch2_OVP, fit.covariate2_OVP,
                         fit.RDmatch3_OVP, fit.covariate3_OVP)

      FE_add = data.frame(matrix(c("Adjustment", "Matching", "Controls", "Matching", "Controls", "Matching", "Controls",
                                   "Sample", "Left", "Left", "Center", "Center", "Right", "Right"),
                                 ncol = length(models_main) + 1, byrow = T))
      attr(FE_add, "position") = c(7, 8)

      regtab = modelsummary(models_main,
                            output = "latex",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "RD results with covariate adjustment. Standard errors in parentheses. Dependent variable: Propensity to vote for OVP \\label{tab:rdtab_ptvOVP_covariates.tex}")
      regtab

  ## Bandwidth selection  ----
    ### Main (Fig A9) ----
      bws = seq(3, 9, 0.5)
      pdat = data.frame()

      for(j in 1:length(bws)){
        est = with(d_combined %>% filter(lr_tri == "left"),
                   rdrobust(ptv_GRUN, running, h = bws[j]))
        est.90 = with(d_combined %>% filter(lr_tri == "left"),
                      rdrobust(ptv_GRUN, running, h = bws[j], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(bw = bws[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "Greens (left voters)"))
      }

      for(j in 1:length(bws)){
        est = with(d_combined %>% filter(lr_tri == "center"),
                   rdrobust(ptv_OVP, running, h = bws[j]))
        est.90 = with(d_combined %>% filter(lr_tri == "center"),
                      rdrobust(ptv_OVP, running, h = bws[j], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(bw = bws[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "ÖVP (centrist voters)"))
      }

      for(j in 1:length(bws)){
        est = with(d_combined %>% filter(lr_tri == "right"),
                   rdrobust(ptv_OVP, running, h = bws[j]))
        est.90 = with(d_combined %>% filter(lr_tri == "right"),
                      rdrobust(ptv_OVP, running, h = bws[j], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(bw = bws[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "ÖVP (right voters)"))
      }

      pdat$selected = F

      pdat$selected[pdat$DV == "Greens (left voters)" & pdat$bw == 3.5] = T
      pdat$selected[pdat$DV == "ÖVP (centrist voters)" & pdat$bw == 3.5] = T
      pdat$selected[pdat$DV == "ÖVP (right voters)" & pdat$bw == 3.5] = T

      p1 = ggplot(pdat, aes(x = bw, y = est, ymin = lower, ymax = upper,
                            group = DV,
                            color = selected))+
        geom_hline(yintercept = 0, linetype = 2) +
        geom_point(size = 1.1,
                   show.legend = F) +
        geom_linerange(linewidth = 0.5,
                       show.legend = F) +
        geom_linerange(inherit.aes = F,
                       mapping = aes(x = bw, y = est, ymin = lower90, ymax = upper90,
                                     group = DV,
                                     color = selected),
                       linewidth = 0.9,
                       show.legend = F) +
        facet_wrap(~DV)  +
        scale_color_manual(values = c("grey50", "black")) +
        scale_y_continuous(limits = c(min(pdat$lower, 0), max(pdat$upper))) +
        theme_bw() +
        theme(panel.grid.minor.x = element_blank(),
              axis.title.x = element_blank(),
              strip.text = element_text(face = "bold", size = 11),
              strip.background = element_blank()) +
        labs(y =  "RD estimate") +
        scale_x_continuous(limits = c(2.8, 9.2), breaks = 3:9)
      p1

      p2 = ggplot(pdat %>%
                    filter(bw %in%1:100) %>%
                    pivot_longer(starts_with("n_"),
                                 names_to = "leftright",
                                 values_to = "n") %>%
                    mutate(leftright = case_when(leftright == "n_pre" ~ "pre treatment",
                                                 leftright == "n_post" ~ "post treatment") %>%
                             factor(., levels = c("pre treatment", "post treatment"))),
                  aes(x = bw, y = n, fill = leftright)) +
        geom_bar(stat = "identity", position = position_dodge(),
                 width = 0.39) +

        scale_fill_manual(values = c("#2E86C1", "#A93226")) +
        theme_bw() +
        theme(legend.title = element_blank(),
              legend.position = "bottom",
              panel.grid.minor.x = element_blank(),
              strip.text = element_blank()) +
        labs(x = "Bandwidth", y =  "Observations") +
        scale_x_continuous(limits = c(2.8, 9.2), breaks = 1:9) +
        scale_y_continuous(expand = expansion(mult = c(0,0.05))) +
        facet_wrap(~DV)
      p2

      ggarrange(p1, p2, ncol = 1,
                heights = c(2, 1))


    ### ATE (Fig A10) ----
      bw = 1:9
      j = 1
      pdat = data.frame()
      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_GRUN ~ treat,
                 data = d_combined %>%
                   filter(lr_tri == "left" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data, pred2$data %>% select(conf.high90 = conf.high,
                                                      conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>%
                           mutate(DV = "Greens (left voters)",
                                  bw = bw[j]))
      }

      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_OVP ~ treat,
                 data = d_combined %>%
                   filter(lr_tri == "center" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data, pred2$data %>% select(conf.high90 = conf.high,
                                                      conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>% mutate(DV = "ÖVP (centrist voters)",
                                         bw = bw[j]))
      }

      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_OVP ~ treat,
                 data = d_combined %>%
                   filter(lr_tri == "right" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data, pred2$data %>% select(conf.high90 = conf.high,
                                                      conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>% mutate(DV = "ÖVP (right voters)",
                                         bw = bw[j]))
      }

      ggplot(pdat, aes(x = bw, y = estimate, ymin = conf.low, ymax = conf.high)) +
        geom_hline(yintercept = 0, linetype = 2) +
        geom_point(size = 2) +
        geom_linerange( linewidth = 0.5,
                        show.legend = F) +
        geom_linerange(inherit.aes = F,
                       mapping = aes(x = bw, y = estimate, ymin = conf.low90, ymax = conf.high90),

                       linewidth = 1,
                       show.legend = F) +
        theme_bw() +
        scale_x_continuous(breaks = pdat$bw) +
        theme(panel.grid.minor.x = element_blank(),
              strip.text = element_text(face = "bold", size = 11),
              strip.background = element_blank()) +
        facet_wrap(~DV) +
        labs(y = "Average treatment effect",
             x = "Bandwidth")

    ### Simple LATE (Fig A11) ----
      bw = 1:9
      j = 1
      pdat = data.frame()
      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_GRUN ~ treat*running,
                 data = d_combined %>%
                   filter(lr_tri == "left" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data[1,], pred2$data[1,] %>% select(conf.high90 = conf.high,
                                                              conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>%
                           mutate(DV = "Greens (left voters)",
                                  bw = bw[j]))
      }

      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_OVP ~ treat*running,
                 data = d_combined %>%
                   filter(lr_tri == "center" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data[1,], pred2$data[1,] %>% select(conf.high90 = conf.high,
                                                              conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>%
                           mutate(DV = "ÖVP (centrist voters)",
                                  bw = bw[j]))
      }

      for(j in 1:length(bw)){
        d_combined$bw = (d_combined$date >= (treatdate - bw[j])) & (d_combined$date <= (treatdate + bw[j] - 1))

        fit = lm(ptv_OVP ~ treat*running,
                 data = d_combined %>%
                   filter(lr_tri == "right" &
                            bw))
        pred = dwplot(fit)
        pred2 = dwplot(fit, ci = 0.9)
        pred = cbind(pred$data[1,], pred2$data[1,] %>% select(conf.high90 = conf.high,
                                                              conf.low90 = conf.low))
        pdat = bind_rows(pdat,
                         pred %>%
                           mutate(DV = "ÖVP (right voters)",
                                  bw = bw[j]))
      }

      ggplot(pdat, aes(x = bw, y = estimate, ymin = conf.low, ymax = conf.high)) +
        geom_hline(yintercept = 0, linetype = 2) +
        geom_point(size = 2) +
        geom_linerange( linewidth = 0.5,
                        show.legend = F) +
        geom_linerange(inherit.aes = F,
                       mapping = aes(x = bw, y = estimate, ymin = conf.low90, ymax = conf.high90),

                       linewidth = 1,
                       show.legend = F) +
        theme_bw() +
        scale_x_continuous(breaks = pdat$bw) +
        theme(panel.grid.minor.x = element_blank(),
              strip.text = element_text(face = "bold", size = 11),
              strip.background = element_blank()) +
        facet_wrap(~DV) +
        labs(y = "Local average treatment effect",
             x = "Bandwidth")

  ## Simple ATE and LATE -----
      summary(fit_main1_regATE <- lm(ptv_GRUN ~ treat, data = d_combined %>% filter(lr_tri == "left")))
      summary(fit_main2_regATE <- lm(ptv_GRUN ~ treat, data = d_combined %>% filter(lr_tri == "center")))
      summary(fit_main3_regATE <- lm(ptv_GRUN ~ treat, data = d_combined %>% filter(lr_tri == "right")))

      summary(fit_main1_regATE_OVP <- lm(ptv_OVP ~ treat, data = d_combined %>% filter(lr_tri == "left")))
      summary(fit_main2_regATE_OVP <- lm(ptv_OVP ~ treat, data = d_combined %>% filter(lr_tri == "center")))
      summary(fit_main3_regATE_OVP <- lm(ptv_OVP ~ treat, data = d_combined %>% filter(lr_tri == "right")))

      summary(fit_main1_regLATE <- lm(ptv_GRUN ~ treat*running, data = d_combined %>% filter(lr_tri == "left")))
      summary(fit_main2_regLATE <- lm(ptv_GRUN ~ treat*running, data = d_combined %>% filter(lr_tri == "center")))
      summary(fit_main3_regLATE <- lm(ptv_GRUN ~ treat*running, data = d_combined %>% filter(lr_tri == "right")))

      summary(fit_main1_regLATE_OVP <- lm(ptv_OVP ~ treat*running, data = d_combined %>% filter(lr_tri == "left")))
      summary(fit_main2_regLATE_OVP <- lm(ptv_OVP ~ treat*running, data = d_combined %>% filter(lr_tri == "center")))
      summary(fit_main3_regLATE_OVP <- lm(ptv_OVP ~ treat*running, data = d_combined %>% filter(lr_tri == "right")))

    ### ATE (Tab A5) ----
      regression_models = list(fit_main1_regATE,
                               fit_main1_regLATE,
                               fit_main2_regATE,
                               fit_main2_regLATE,
                               fit_main3_regATE,
                               fit_main3_regLATE)

      FE_add = data.frame(matrix(c("Estimator", "ATE", "LATE", "ATE", "LATE", "ATE", "LATE",
                                   "Sample", "Left", "Left", "Center", "Center", "Right", "Right"),
                                 ncol = length(regression_models) + 1, byrow = T))
      attr(FE_add, "position") = c(9, 10)

      regtab = modelsummary(regression_models,
                            output = "latex",
                            add_rows = FE_add,
                            coef_map = cmap,
                            escape = F,
                            stars = T,
                            gof_omit = "AIC|BIC|Log|F|RMSE",
                            title = "Linear regression estimates. Standard errors in parentheses. Dependent variable: Propensity to vote for Greens \\label{tab:regtab_ptvGRUN}")
      regtab

    ### LATE (Tab A6) ----
      regression_models = list(fit_main1_regATE_OVP,
                               fit_main1_regLATE_OVP,
                               fit_main2_regATE_OVP,
                               fit_main2_regLATE_OVP,
                               fit_main3_regATE_OVP,
                               fit_main3_regLATE_OVP)

      FE_add = data.frame(matrix(c("Estimator", "ATE", "LATE", "ATE", "LATE", "ATE", "LATE",
                                   "Sample", "Left", "Left", "Center", "Center", "Right", "Right"),
                                 ncol = length(regression_models) + 1, byrow = T))
      attr(FE_add, "position") = c(9, 10)

      regtab = modelsummary(regression_models,
                            output = "latex",
                            add_rows = FE_add,
                            coef_map = cmap,
                            escape = F,
                            stars = T,
                            gof_omit = "AIC|BIC|Log|F|RMSE",
                            title = "Linear regression estimates. Standard errors in parentheses. Dependent variable: Propensity to vote for ÖVP \\label{tab:regtab_ptvOVP}")
      regtab

  ## Excluding Sunday ----
    ### RDD plot (Fig A12) ----
      d_combined_alt = d_combined %>%
        mutate(running = as.numeric(date - as.Date("2024-06-16"))) %>%
        filter(running != 0) %>%
        group_by(lr_tri) %>%
        mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T),
               ptv_OVP = ptv_OVP - mean(ptv_OVP,na.rm = T))

      grid.arrange(RDDplot(d_combined_alt, "ptv_GRUN", "running", "treat", "lr_tri_label") +
                     scale_y_continuous(limits = c(-2.4, 1.5)) +
                     labs(y = "Propensity to vote for Greens (centered)",
                          x = " "),
                   RDDplot(d_combined_alt, "ptv_OVP", "running", "treat", "lr_tri_label") +
                     scale_y_continuous(limits = c(-2.4, 1.5))+
                     labs(y = "Propensity to vote for ÖVP (centered)",
                          x = "Days since Press Conference"))

    ### RDD table (Tab A7) ----
      summary(fit1 <- with(d_combined_alt %>% filter(lr_tri == "left"),
                           rdrobust(ptv_GRUN, running, all = T)))
      summary(fit2 <- with(d_combined_alt %>% filter(lr_tri == "center"),
                           rdrobust(ptv_GRUN, running, all = T)))
      summary(fit3 <- with(d_combined_alt %>% filter(lr_tri == "right"),
                           rdrobust(ptv_GRUN, running, all = T)))

      summary(fit1_OVP <- with(d_combined_alt %>% filter(lr_tri == "left"),
                               rdrobust(ptv_OVP, running, all = T)))
      summary(fit2_OVP <- with(d_combined_alt %>% filter(lr_tri == "center"),
                               rdrobust(ptv_OVP, running, all = T)))
      summary(fit3_OVP <- with(d_combined_alt %>% filter(lr_tri == "right"),
                               rdrobust(ptv_OVP, running, all = T)))
      models_main = list(fit1, fit2, fit3,
                         fit1_OVP, fit2_OVP, fit3_OVP)

      FE_add = data.frame(matrix(c("Sample", "Left", "Center", "Right",
                                   "Left", "Center", "Right"),
                                 ncol = length(models_main) + 1, byrow = T))
      attr(FE_add, "position") = c(7, 8)

      regtab = modelsummary(models_main,
                            output = "latex",
                            add_rows = FE_add,
                            escape=FALSE,
                            stars = T,
                            title = "RD results excluding 16 June. Standard errors in parentheses. Dependent variables: Propensity to vote for Greens (models 1-3) and ÖVP (models 4-6). \\label{tab:rdtab_ptvGRUNOVP_excludeSunday}")# c('*' = .05, '**' = .01, "***" = .001))
      regtab

  ## Placebo ----
    ### Pre-treatment dates (Fig A13) ----
      d_combined_placebo = d_combined
      min(d_combined$date)
      placebo_dates = seq(as.Date("2024-06-14"), as.Date("2024-06-17"), 1)
      pdat = data.frame()

      for(j in 1:length(placebo_dates)){
        d_combined_placebo = d_combined_placebo %>%
          mutate(treat = date >= placebo_dates[j],
                 running = as.numeric(date - placebo_dates[j]) + 0.5)
        est = with(d_combined_placebo %>% filter(lr_tri == "left"),
                   rdrobust(ptv_GRUN, running, h = fit1.main$bws[1,1]))
        est.90 = with(d_combined_placebo %>% filter(lr_tri == "left"),
                      rdrobust(ptv_GRUN, running, h = fit1.main$bws[1,1], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(placebo_date = placebo_dates[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "Greens (left voters)"))
      }

      for(j in 1:length(placebo_dates)){
        d_combined_placebo = d_combined_placebo %>%
          mutate(treat = date >= placebo_dates[j],
                 running = as.numeric(date - placebo_dates[j]) + 0.5)
        est = with(d_combined_placebo %>% filter(lr_tri == "center"),
                   rdrobust(ptv_OVP, running, h = fit2.main_OVP$bws[1,1]))
        est.90 = with(d_combined_placebo %>% filter(lr_tri == "center"),
                      rdrobust(ptv_OVP, running, h = fit2.main_OVP$bws[1,1], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(placebo_date = placebo_dates[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "ÖVP (centrist voters)"))
      }

      for(j in 1:length(placebo_dates)){
        d_combined_placebo = d_combined_placebo %>%
          mutate(treat = date >= placebo_dates[j],
                 running = as.numeric(date - placebo_dates[j]) + 0.5)
        est = with(d_combined_placebo %>% filter(lr_tri == "right"),
                   rdrobust(ptv_OVP, running, h = fit3.main_OVP$bws[1,1]))
        est.90 = with(d_combined_placebo %>% filter(lr_tri == "right"),
                      rdrobust(ptv_OVP, running, h = fit3.main_OVP$bws[1,1], level = 90))

        pdat = bind_rows(pdat,
                         data.frame(placebo_date = placebo_dates[j],
                                    est = est$Estimate[1],
                                    lower = est$ci[1,1],
                                    upper = est$ci[1,2],
                                    lower90 = est.90$ci[1,1],
                                    upper90 = est.90$ci[1,2],
                                    est.robust = est$Estimate[2],
                                    lower.robust = est$ci[3,1],
                                    upper.robust = est$ci[3,2],
                                    n_pre = est$N_h[1],
                                    n_post = est$N_h[2],
                                    DV = "ÖVP (right voters)"))
      }

      pdat$selected = (pdat$placebo_date == treatdate)


      ggplot(pdat, aes(x = placebo_date, y = est, ymin = lower, ymax = upper, color = selected))+
        geom_hline(yintercept = 0, linetype = 2) +
        geom_pointrange(size = 1,
                        linewidth = 1,
                        show.legend = F) +
        geom_linerange(inherit.aes = F,
                       mapping = aes(x = placebo_date, y = est, ymin = lower90, ymax = upper90,
                                     color = selected),
                       linewidth = 2,
                       show.legend = F) +
        scale_color_manual(values = c("grey50", "black")) +
        scale_y_continuous(limits = c(min(pdat$lower, 0), max(pdat$upper))) +
        theme_bw() +
        theme(panel.grid.minor.x = element_blank(),
              axis.title.x = element_blank(),
              strip.text = element_text(face = "bold", size = 11),
              strip.background = element_blank(),
              axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1.1)) +
        scale_x_date(breaks = pdat$placebo_date[pdat$selected] + seq(-4, 4, 1),
                     date_minor_breaks = "day",
                     date_labels = "%d %B") +
        labs(y =  "RD estimate",
             x = "Placebo date") +
        facet_wrap(~ DV)

    ### All AUTNES waves (Fig A14) ----
      ests = data.frame()
      d_autnes_all = d_autnes_all %>% group_by(wave) %>%
        mutate(treat_date = median(date[date != names(table(date))[which.max(table(date))]], na.rm = T),
               treat = date >= treat_date,
               running = as.numeric(date - treat_date) + 0.5)

      waves = unique(d_autnes_all$wave)
      for(j in 2:length(waves)){
        d = d_autnes_all %>% filter(wave == waves[j] & lr_tri == "left")
        est = try(rdrobust(d$ptv_GRUN, d$running, h = 4))
        est90 = try(rdrobust(d$ptv_GRUN, d$running, h = 4, level = 90))

        if(!inherits(est, "try-error"))
          ests = bind_rows(ests,
                           data.frame(wave = as.character(waves[j]),
                                      est = est$Estimate[1],
                                      lower = est$ci[1,1],
                                      upper = est$ci[1,2],
                                      lower90 = est90$ci[1,1],
                                      upper90 = est90$ci[1,2],
                                      est.robust = est$Estimate[2],
                                      lower.robust = est$ci[3,1],
                                      upper.robust = est$ci[3,2],
                                      n_pre = est$N_h[1],
                                      n_post = est$N_h[2]))

      }

      # add main finding
      summary(fit1.main <- with(d_combined %>% filter(lr_tri == "left"),
                                rdrobust(ptv_GRUN, running, all = T)))
      summary(fit1.main90 <- with(d_combined %>% filter(lr_tri == "left"),
                                  rdrobust(ptv_GRUN, running, all = T, level = 90)))

      ests$actual = F
      ests = bind_rows(ests,
                       data.frame(wave = "Main finding",
                                  actual = T,
                                  est = fit1.main$Estimate[1],
                                  lower = fit1.main$ci[1,1],
                                  upper = fit1.main$ci[1,2],
                                  lower90 = fit1.main90$ci[1,1],
                                  upper90 = fit1.main90$ci[1,2],
                                  est.robust = fit1.main$Estimate[2],
                                  lower.robust = fit1.main$ci[3,1],
                                  upper.robust = fit1.main$ci[3,2],
                                  n_pre = fit1.main$N_h[1],
                                  n_post = fit1.main$N_h[2]))
      ests = ests %>%
        arrange(est) %>%
        mutate(order = 1:n(),
               wave = factor(wave, levels = unique(wave)))


      ggplot(ests, aes(x = est, xmin = lower, xmax = upper, y = wave, color = actual)) +
        geom_vline(xintercept = 0, linetype = 2) +
        geom_pointrange(position = position_dodge(width = 0.6),
                        size = 1,
                        linewidth = 1,
                        show.legend = F) +
        geom_linerange(inherit.aes = F,
                       mapping = aes(x = est, xmin = lower90, xmax = upper90,
                                     y = wave, color = actual),
                       linewidth = 2,

                       position = position_dodge(width = 0.6),
                       show.legend = F) +
        #facet_wrap(~ actual, ncol = 1, drop = T, scales = "free_y") +
        theme_bw() +
        scale_color_manual(values = c("grey50", "black")) +
        theme(legend.title = element_blank(),
              legend.position = "bottom",
              legend.key.width = unit(2, "cm"),
              axis.title.y = element_blank()) +
        labs(x = "RD estimate")

  ## Manipulation checks ----
    ### MIP (Fig A15) ----
      d_combined_plot = d_combined %>%
        filter(!is.na(lr_tri))
      p1 = RDDplot(d_combined_plot, "mip_climateenvironment", "running", "treat") +
        labs(y = "Climate change / environment\namong most important issues")
      p2 = RDDplot(d_combined_plot, "mip_immigration", "running", "treat") +
        labs(y = "Immigration\namong most important issues")
      ggarrange(p1 + coord_cartesian(ylim = c(0.125, 0.375)),
                p2 + coord_cartesian(ylim = c(0.125, 0.375)),
                ncol = 2)

    ### Fear (Fig A16) ----
      d_combined_plot = d_combined %>%
        filter(!is.na(lr_tri)) %>%
        mutate(fear_climatechange = fear_climatechange - mean(fear_climatechange, na.rm = T),
               fear_immigration = fear_immigration - mean(fear_immigration, na.rm = T))
      p1 = RDDplot(d_combined_plot, "fear_climatechange", "running", "treat") +
        labs(y = "Fear of climate change (centered)")
      p2 = RDDplot(d_combined_plot, "fear_immigration", "running", "treat") +
        labs(y = "Fear of immigration (centered)")
      ggarrange(p1 + coord_cartesian(ylim = c(-0.7, 0.7)),
                p2 + coord_cartesian(ylim = c(-0.7, 0.7)),
                ncol = 2)

  ## Non-response (Fig A17) ----
    d_combined_plot = d_combined %>%
      filter(!is.na(lr_tri)) %>%
      mutate(nonresponse_ptvGRUN = as.numeric(is.na(ptv_GRUN)),
             nonresponse_ptvOVP = as.numeric(is.na(ptv_OVP)))


    grid.arrange(RDDplot(d_combined_plot, "nonresponse_ptvGRUN", "running", "treat", "lr_tri_label") +
                   #scale_y_continuous(limits = c(-2.4, 1.5)) +
                   labs(y = "Nonresponse PTV for Greens",
                        x = " "),
                 RDDplot(d_combined_plot, "nonresponse_ptvOVP", "running", "treat", "lr_tri_label") +
                   #scale_y_continuous(limits = c(-2.4, 1.5))+
                   labs(y = "Nonresponse PTV for ÖVP"))

  ## Seoarate surveys ----
    ### Greens (Fig A18) ----
      d_combined_plot = d_combined %>%
        filter(survey == "ees") %>%
        group_by(lr_tri) %>%
        mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T))
      p1 = RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "lr_tri_label") +
        theme(plot.title = element_text(face = "bold")) +
        labs(y = "Propensity to vote for Greens (centered)",
             x = "Days since Council vote",
             title = "EES")

      d_combined_plot = d_combined %>%
        filter(survey == "autnes") %>%
        group_by(lr_tri) %>%
        mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T))
      p2 = RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "lr_tri_label") +
        scale_y_continuous(breaks = c(-2.0, -1.0, 0, 1.0, 2.0)) +
        theme(plot.title = element_text(face = "bold")) +
        labs(y = "Propensity to vote for Greens (centered)",
             x = "Days since Council vote",
             title = "AUTNES")
      p3 = ggarrange(p1 + rremove("ylab"),
                     p2 + rremove("ylab"),
                     ncol = 1)

      annotate_figure(p3,
                      left = textGrob("Propensity to vote for Greens (centered)", rot = 90, vjust = 1, gp = gpar(cex = 1)))

    ### OVP (Fig A19) ----
      d_combined_plot = d_combined %>%
        filter(survey == "ees") %>%
        group_by(lr_tri) %>%
        mutate(ptv_OVP = ptv_OVP - mean(ptv_OVP,na.rm = T))
      p1 = RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "lr_tri_label") +
        scale_y_continuous(breaks = c(-4.0, -2.0, 0, 2.0, 4.0)) +
        theme(plot.title = element_text(face = "bold")) +
        labs(y = "Propensity to vote for ÖVP (centered)",
             x = "Days since Council vote",
             title = "EES")

      d_combined_plot = d_combined %>%
        filter(survey == "autnes") %>%
        group_by(lr_tri) %>%
        mutate(ptv_OVP = ptv_OVP - mean(ptv_OVP,na.rm = T))
      p2 = RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "lr_tri_label") +
        theme(plot.title = element_text(face = "bold")) +
        labs(y = "Propensity to vote for ÖVP (centered)",
             x = "Days since Council vote",
             title = "AUTNES")
      p3 = ggarrange(p1 + rremove("ylab"),
                     p2 + rremove("ylab"),
                     ncol = 1)

      annotate_figure(p3,
                      left = textGrob("Propensity to vote for ÖVP (centered)", rot = 90, vjust = 1, gp = gpar(cex = 1)))


  ## Position relative to compromise (Fig A20) ----
    d_combined_plot = d_combined %>%
      group_by(leftofcompromise_label) %>%
      mutate(ptv_GRUN = ptv_GRUN - mean(ptv_GRUN,na.rm = T),
             ptv_OVP = ptv_OVP - mean(ptv_OVP, na.rm = T))

    grid.arrange(RDDplot(d_combined_plot, "ptv_GRUN", "running", "treat", "leftofcompromise_label") +
                   scale_y_continuous(limits = c(-1.5, 1.25)) +
                   labs(y = "Propensity to vote for Greens (centered)",
                        x = " "),
                 RDDplot(d_combined_plot, "ptv_OVP", "running", "treat", "leftofcompromise_label") +
                   scale_y_continuous(limits = c(-1.5, 1.25))+
                   labs(y = "Propensity to vote for ÖVP (centered)"))
